home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 3 / BBS in a box - Trilogy III.iso / Files / Prog / U-Z / VideoToolBox Folder / Utilities / Quick3 / PsychometricFit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-02-23  |  6.2 KB  |  171 lines  |  [TEXT/KAHL]
  1. /* 
  2. PsychometricFit.c
  3. Copyright 1990 (c) Denis G. Pelli
  4. A general-purpose function that does a maximum likelihood fit of any psychometric
  5. function to psychometric data. The returned value is the level of significance
  6. at which the fit can be rejected. The degreesOfFreedom may be zero, in which case 
  7. no parameters
  8. will be adjusted, but you'll get the log likelihood and significance of the supplied
  9. parameter values.
  10.  
  11. The psychometric function (which you supply as an argument) takes two arguments:
  12. a contrast and a pointer to a paramRecord.
  13. The function Weibull() is provided in Weibull.c. Others may be added as desired.
  14. It is assumed that the first parameter is the log of threshold. No assumptions are
  15. made about the other parameters, except that any that are to be iteratively fit are
  16. assumed to be of type "double".
  17.  
  18. Quick3 is a stand-alone program that uses PsychometricFit() to do the real work.
  19. I suggest you read the introductory comments at the beginning of Quick3.c
  20.  
  21. HISTORY:
  22. 4/7/90        dgp wrote it
  23. 10/29/90    dgp    tidied up the comments
  24. 11/17/92    dgp "
  25. 1/25/93 dgp removed obsolete support for THINK C 4.
  26.  
  27. SOURCES:
  28. Quick3.h
  29. LogLikelihood.c
  30. MonotonicFit.c
  31. PsychometricFit.c
  32. SortAndMergeContrasts.c
  33. Weibull.c
  34. #From Denis Pelli's VideoToolbox:
  35. VideoToolbox.h
  36. Binomial.c
  37. ChiSquare.c
  38. MyFgets.c
  39. Normal.c
  40. SetFileInfo.c        # Used only on the Macintosh
  41. #From Numerical Recipes in C:
  42. nr.h
  43. NRUTIL.h
  44. BRENT.C
  45. F1DIM.C
  46. LINMIN.C
  47. MNBRAK.C
  48. NRUTIL.C
  49. POWELL.C
  50.  
  51. LIMITATIONS
  52.  
  53. This program uses routines from Numerical Recipes in C. They're copyrighted, 
  54. so I can't distribute them. You can order the software:
  55.     Numerical Recipes C Diskette for Macintosh $29.95
  56. and book:
  57.     Numerical Recipes in C: The Art of Scientific Computing $44.50
  58. from:
  59.     Cambridge University Press
  60.     Order Department
  61.     510 North Avenue
  62.     New Rochelle, NY 10801
  63.  
  64. Note that I have made several changes to the Numerical Recipes in C routines: 
  65. 1.In every file I changed "float" to "FLOAT", and #included nr.h. I inserted the 
  66. statement "typedef double FLOAT;" in the file nr.h. This is because the 
  67. Macintosh computes doubles much faster than floats. If you'd rather run
  68. slowly than modify your Numerical Recipes in C files, then you will need to insert 
  69. "typedef float FLOAT;"
  70. in CalibrateLuminance.c & PsychometricFit.c in order to compile those files.
  71. The rest of the VideoToolbox doesn't care.
  72. 2. See the Read me! document for other changes and bug fixs.
  73. */
  74. #include "VideoToolbox.h"
  75. #include "Quick3.h"
  76. #include <nr.h>
  77.  
  78. #define TOLERANCE 0.001    /* fractional tolerance of log likelihood. Not critical. */
  79.  
  80. /*
  81. I hate global variables because they hide the flow of information. However,
  82. some sort of cludge is necessary to pass the extra arguments to Error(), bypassing
  83. the Numerical Recipes routines that call it, since they only pass the 
  84. parameters they know about. These static declations at least restrict the scope of
  85. these "globals" to this file.
  86. */
  87. static double Error(FLOAT *p);
  88. static dataRecord *myDataPtr;                /* for Error() */
  89. static PsychometricFunctionPtr MyPsychFun;    /* for Error() */
  90. static paramRecord myParams;                /* for Error() */
  91. static int myDegreesOfFreedom;                /* for Error() */
  92. static int iter;                            /* for Error() */
  93.  
  94. double PsychometricFit(paramRecord *paramPtr,PsychometricFunctionPtr PsychFun
  95.     ,dataRecord *dataPtr,double *logLikelihoodPtr,int degreesOfFreedom
  96.     ,double *chiSquarePtr,int *chiSquareDFPtr)
  97. {
  98.     int i,j;
  99.     FLOAT *p,**direction,ftol,fret;
  100.     dataRecord monotonicData;
  101.     double monotonicLL;
  102.     int monotonicDF;
  103.     double P;
  104.     
  105.     myDataPtr=dataPtr;        /* copy these for use by Error() */
  106.     MyPsychFun=PsychFun;
  107.     myParams=*paramPtr;
  108.     myDegreesOfFreedom=degreesOfFreedom;
  109.     
  110.     p=vector(1,degreesOfFreedom);
  111.     direction=matrix(1,degreesOfFreedom,1,degreesOfFreedom);    /* initial set of directions */
  112.     if(p==NULL || direction == NULL) {
  113.         PrintfExit("PsychometricFit: not enough room for arrays.\007\n");
  114.     }
  115.  
  116.     for(i=1;i<=degreesOfFreedom;i++) p[i]=((double *)paramPtr)[i-1];
  117.     for(i=1;i<=degreesOfFreedom;i++)for(j=1;j<=degreesOfFreedom;j++)direction[i][j]=0.0;
  118.     for(i=1;i<=degreesOfFreedom;i++)direction[i][i]=0.03;    /* initial step size */
  119.     ftol=TOLERANCE;    /* fractional tolerance on Error value when done */
  120.     iter=0;
  121.     
  122.     /* do it. The psychomtric function is passed to Error by the global MyPsychFun */
  123.     if(degreesOfFreedom==0) fret = Error(p);
  124.     else powell(p,direction,degreesOfFreedom,ftol,&iter,&fret,&Error);
  125.  
  126.     for(i=1;i<=degreesOfFreedom;i++) ((double *)paramPtr)[i-1]=p[i];
  127.     free_matrix(direction,1,degreesOfFreedom,1,degreesOfFreedom);
  128.     free_vector(p,1,degreesOfFreedom);
  129.  
  130.     *logLikelihoodPtr = -fret;
  131.     
  132.     /* Now compute the degree of significance at which the fit can be rejected */
  133.     monotonicData= *dataPtr;
  134.     MonotonicFit(&monotonicData,&monotonicLL,&monotonicDF);    /* overwrites data with fit */
  135.     *chiSquarePtr= -2.0*(*logLikelihoodPtr-monotonicLL);    /* -2 log likelihood ratio of hypotheses */
  136.     *chiSquareDFPtr=monotonicDF-degreesOfFreedom;            /* difference in degrees of freedom */
  137.     P=PChiSquare(*chiSquarePtr,*chiSquareDFPtr);            /* significance */
  138.     return P;
  139. }
  140.  
  141. /* There is a subtlety here. I thought that I could use Powell with the whole
  142. paramRecord, yet ask Powell to only twiddle the first few parameters, figuring that
  143. even when I was asking Powell to fit only the first few parameters the other
  144. parameters would still be there in the array when the pointer to the array was
  145. passed to Error(). Alas, Powell() and its subroutines make COPIES of the array,
  146. and naturally fail to copy the non-twiddled parameters, since they don't know about
  147. them. The solution is for Error() to have its own complete copy of the paramRecord.
  148. Each time Error() is called it updates the twiddled parameters before calling
  149. LogLikelihood(), which calls the psychometric function (*MyPsychFun)().
  150. */
  151.  
  152. double Error(FLOAT *p)
  153. {
  154.     double error;
  155.     int i;
  156.     static int lastIter=0;
  157.     
  158.     for(i=1;i<=myDegreesOfFreedom;i++) ((double *)&myParams)[i-1]=p[i];
  159.     
  160.     error = -LogLikelihood(myDataPtr,&myParams,MyPsychFun);
  161.     
  162.     /* Diagnostic printout for difficult cases */
  163.     if(iter>0 && iter%50 == 0 && iter!=lastIter){
  164.         printf("Error(): Warning, %d iterations:\n",iter);
  165.         printf("logAlpha %5.2f, beta %5.1f, gamma %5.2f, delta %6.3f -log likelihood %9.0g\n",myParams.logAlpha,myParams.beta,myParams.gamma,myParams.delta,error);
  166.         lastIter=iter;
  167.     }
  168.     return error;
  169. }
  170.  
  171.